from copy      import deepcopy
from re        import search as re_search
from tempfile  import gettempdir
from itertools import product as itertools_product
import collections
import numpy 
import sys

# --------------------------------------------------------------------
#   A dictionary of useful constants.
#
#   Whilst the dictionary may be modified directly, it is safer to
#   retrieve and set the values with a function where one is
#   provided. This is due to interdependencies between some values.
#
#   :Keys:
#
#        ATOL : float
#	    The value of absolute tolerance for testing numerically
#	    tolerant equality. Retrieved and set with :func:`.ATOL`.
#
#        CHUNKSIZE : int
#	    The chunk size (in bytes) for data storage and
#	    processing. Retrieved and set with :func:`.CHUNKSIZE`.
#
#        FM_THRESHOLD : float
#	    The minimum amount of memory (in kibibytes) to be kept
#	    free for temporary work space. Retrieved and set with
#	    :func:`.FM_THRESHOLD`.
#
#        MINNCFM : int
#	    The number of chunk sizes to be kept free for temporary
#	    work space.
#
#        RTOL : float
#	    The value of relative tolerance for testing numerically
#	    tolerant equality. Retrieved and set with :func:`.RTOL`.
#
#        TEMPDIR : str
#	    The location to store temporary files. By default, use the
#	    default directory used by the :mod:`tempfile` module.
# --------------------------------------------------------------------
CONSTANTS = {'RTOL'   : sys.float_info.epsilon,
             'ATOL'   : sys.float_info.epsilon,
             'TEMPDIR': gettempdir(),
             }

_open_files = []

def ATOL(*atol):
    '''

Return or set the value of absolute tolerance for testing numerically
tolerant equality.

:Parameters:

    atol : int, optional
        The new value of absolute tolerance.

:Returns:

    out : float or None    
        If `atol` was not set return the existing value of absolute
        tolerance, otherwise return `None`.

**Examples**

>>> cf.ATOL()
1e-08
>>> cf.ATOL(1e-10)
>>> cf.ATOL()
1e-10

'''
    if atol:
        CONSTANTS['ATOL'] = atol[0]
    else:
        return CONSTANTS['ATOL']
#--- End: def


# ====================================================================
#
# CfList object
#
# ====================================================================

class CfList(collections.MutableSequence):
    '''

A list-like object with attributes.

'''
    def __init__(self, sequence=()):
        '''

**Initialization**

:Parameters:

    sequence : iterable, optional
         Define a new list with these elements.

'''
        self.__dict__['_list'] = []

        try:
            self.__dict__['_list'].extend(sequence)
        except TypeError:
            self.__dict__['_list'].append(sequence)
    #--- End: def

    def __deepcopy__(self, memo):
        '''
Used if copy.deepcopy is called on the variable.

''' 
        return self.copy()
    #--- End: def

    def __repr__(self):
        '''
x.__repr__() <==> repr(x)

'''
        return repr(self._list)
    #--- End: def

    def __str__(self):
        '''
x.__str__() <==> str(x)

'''
        return str(self._list)
    #--- End: def

    def __contains__(self, item):
        '''
x.__contains__(y) <==> y in x

Uses numerically tolerant equality where appropriate.

'''
        for x in self._list:
            if equals(item, x, traceback=False):
                return True
        #--- End: for

        return False
    #--- End: def
    
    def __len__(self):
        '''
x.__len__() <==> len(x)

'''
        return len(self._list)
    #--- End: def

    def __delitem__(self, index):
        '''
x.__delitem__(index) <==> del x[index]

'''
        del self._list[index]
    #--- End: def

    def __getitem__(self, index):
        '''
x.__getitem__(index) <==> x[index]

'''
        out = self._list[index]
        try:
            int(index)
        except TypeError:
            # index is a slice instance, so return a CfList instance
            out = type(self)(out)

        return out
    #--- End: def

    def __setitem__(self, index, value):        
        '''
x.__setitem__(index, value) <==> x[index]=value

'''
        try:
            self._list[index] = value._list
        except AttributeError:
            self._list[index] = value
    #--- End: def

    def __eq__(self, other):
        '''
x.__eq__(y) <==> x==y <==> x.equals(y)

'''
        return self.equals(other)
    #--- End: def

    def __ne__(self, other):
        '''
x.__ne__(y) <==> x!=y <==> not x.equals(y)

'''
        return not self.__eq__(other)
    #--- End: def

    def __add__(self, other):
        '''
x.__add__(y) <==> x+y

'''
        new = self.copy()
        new.extend(other)
        return new
    #--- End: def

    def __mul__(self, other):
        '''
x.__mul__(n) <==> x*n

'''
        return type(self)(self._list * other)
    #--- End: def

    def __rmul__(self, other):
        '''
x.__rmul__(n) <==> x*n

'''
        return self.__mul__(other)
    #--- End: def

    def __iadd__(self, other):
        '''
x.__iadd__(y) <==> x+=y

'''
        self.extend(other)
        return self
    #--- End: def

    def __imul__(self, other):
        '''
x.__imul__(n) <==> x*=n

'''
        self._list = self._list * other
        return self
    #--- End: def

    def count(self, value):
        '''

Return the number of occurrences of a given value.

Uses numerically tolerant equality where appropriate.

:Parameters:

    value :
        The value to count.

:Returns:

    out : int
        The number of occurrences of `value`.

**Examples**

>>> s
[1, 2, 3, 2, 4, 2]
>>> s.count(1)
1
>>> s.count(2)
3

'''
        return sum(1
                   for x in self._list 
                   if equals(value, x, traceback=False))
    #--- End def

    def index(self, value, start=0, stop=None):
        '''

Return the first index of a given value.

Uses numerically tolerant equality where appropriate.

:Parameters:

    value :
        The value to look for in the list.

    start : int, optional
        Start looking from this index. By default, look from the
        beginning of the list.

    stop : int, optional
        Stop looking before this index. By default, look up to the end
        of the list.

:Returns:

    out : int

:Raises:

    ValueError :
        If the given value is not in the list.

**Examples**

>>> s
[1, 2, 3, 2, 4, 2]
>>> s.index(1)
1
>>> s.index(2, 2)
3
>>> s.index(2, start=2, stop=5)
3
>>> s.index(6)
ValueError: CfList doesn't contain: 6

'''      
        if start < 0:
            start = len(self) + start

        if stop is None:
            stop = len(self)
        elif stop < 0:
            stop = len(self) + stop

        for i, x in enumerate(self[start:stop]):
            if equals(value, x, traceback=False):
               return i + start
        #--- End: for

        raise ValueError("%s doesn't contain: %s" % 
                         (self.__class__.__name__, repr(value)))
    #--- End: def

    def insert(self, index, item):
        '''

Insert an object before the given index in place.

:Parameters:

    index : int

    item :

:Returns:

    None

**Examples**

>>> s
[1, 2, 3]
>>> s.insert(1, 'A')
>>> s
[1, 'A', 2, 3]
>>> s.insert(100, 'B')
[1, 'A', 2, 3, 'B']
>>> s.insert(100, 'B')
[1, 'A', 2, 3, 'B']
>>> s.insert(-2, 'C')
[1, 'A', 2, 'C', 3, 'B']
>>> s.insert(-100, 'D')
['D', 1, 'A', 2, 'C', 3, 'B']

'''
        if hasattr(item, '_list'):
            self._list.insert(index, item._list)
        else:
            self._list.insert(index, item)
    #--- End: def

    def copy(self):
        '''

Return a deep copy.

Equivalent to ``copy.deepcopy(s)``.

:Returns:

    out : 
        The deep copy.

**Examples**

>>> s.copy()

'''
        new = type(self)()

        for attr, value in self.__dict__.iteritems():
            setattr(new, attr, deepcopy(value))

        return new
    #--- End: def

    def equals(self, other, rtol=None, atol=None, traceback=False):
        '''

True if two instances are equal, False otherwise.

Two instances are equal if their attributes are equal and their
elements are equal pair-wise.

:Parameters:

    other : 
        The object to compare for equality.

    atol : float, optional
        The absolute tolerance for all numerical comparisons, By
        default the value returned by the `ATOL` function is used.

    rtol : float, optional
        The relative tolerance for all numerical comparisons, By
        default the value returned by the `RTOL` function is used.

    traceback : bool, optional
        If True then print a traceback highlighting where the two
        instances differ.

:Returns: 

    out : bool
        Whether or not the two instances are equal.

'''
        if self is other:
            return True
         
        # Check that each instance is the same type
        if self.__class__ != other.__class__:
            if traceback:
                print("%s: Different types: %s, %s" %
                      (self.__class__.__name__,
                       self.__class__.__name__,
                       other.__class__.__name__))
            return False
        #--- End: if

        # Check that the lists have the same number of elements
        if len(self) != len(other): 
            if traceback:
                print("%s: Different attributes: %s, %s" %
                      (self.__class__.__name__,
                       attrs.symmetric_difference(other.__dict__)))
            return False
        #--- End: if

        # Check the attributes
        attrs = set(self.__dict__)
        if attrs != set(other.__dict__):
            if traceback:
                print("%s: Different attributes: %s, %s" %
                      (self.__class__.__name__,
                       attrs.symmetric_difference(other.__dict__)))
            return False
        #--- End: if

        if rtol is None:
            rtol = RTOL()
        if atol is None:
            atol = ATOL()

        for attr, value in attrs - set(('_list',)):
            x = getattr(self, attr)
            y = getattr(other, attr)
            if not equals(x, y, rtol=rtol, atol=atol, traceback=traceback):
                if traceback:
                    print("%s: Different '%s': %s, %s" %
                          (self.__class__.__name__, attr, x, y))
                return False
        #--- End: for

        # Check the element values
        for x, y in zip(self._list, other._list):
            if not equals(x, y, rtol=rtol, atol=atol, traceback=traceback):
                if traceback:
                    print("%s: Different elements: %s, %s" %
                          (self.__class__.__name__, repr(x), repr(y)))
                return False
        #--- End: for

        return True
    #--- End: def
    
#--- End: class


# ====================================================================
#
# CfDict object
#
# ====================================================================

class CfDict(collections.MutableMapping):
    '''

A dictionary-like object with attributes.

'''
    def __init__(self, *args, **kwargs):
        '''

**Initialization**

:Parameters:

    args, kwargs :
        Keys and values are initialized exactly as for a built-in
        dict.

'''

        self._dict = dict(*args, **kwargs)
    #--- End: def

    def __deepcopy__(self, memo):
        '''
Used if copy.deepcopy is called on the variable.

''' 
        return self.copy()
    #--- End: def

    def __repr__(self):
        '''
x.__repr__() <==> repr(x)

'''
        return repr(self._dict)
    #--- End: def

    def __str__(self):
        '''
x.__str__() <==> str(x)

'''
        return str(self._dict)
    #--- End: def

    def __getitem__(self, key):
        '''
x.__getitem__(key) <==> x[key]

'''     
        return self._dict[key]
    #--- End: def

    def __setitem__(self, key, value):
        '''
x.__setitem__(key, value) <==> x[key]=value

'''
        self._dict[key] = value
    #--- End: def

    def __delitem__(self, key):
        '''
x.__delitem__(key) <==> del x[key]

'''
        del self._dict[key]
    #--- End: def

    def __iter__(self):
        '''
x.__iter__() <==> iter(x)

'''
        return iter(self._dict)
    #--- End: def

    def __len__(self):
        '''
x.__len__() <==> len(x)

'''
        return len(self._dict)
    #--- End: def

#    def __contains__(self, item):
#        '''
#Test for set membership using numerically tolerant equality.
#'''
#        for s in self:
#            if equals(item, s):
#                return True
#
#        return False
#    #--- End: def

    def __eq__(self, other):
        '''
x.__eq__(y) <==> x==y <==> x.equals(y)

'''
        return self.equals(other)
    #--- End: def

    def __ne__(self, other):
        '''
x.__ne__(y) <==> x!=y <==> not x.equals(y)

'''
        return not self.__eq__(other)
    #--- End: def

    def copy(self):
        '''

Return a deep copy.

Equivalent to ``copy.deepcopy(d)``.

:Returns:

    out : 
        The deep copy.

**Examples**

>>> d.copy()

'''
        new = type(self)()

        for attr, value in self.__dict__.iteritems():
            setattr(new, attr, deepcopy(value))

        return new
    #--- End: def

    def equals(self, other, rtol=None, atol=None, traceback=False):
        '''

True if two instances are logically equal, False otherwise.

:Parameters:

    other : 
        The object to compare for equality.

    atol : float, optional
        The absolute tolerance for all numerical comparisons, By
        default the value returned by the `ATOL` function is used.

    rtol : float, optional
        The relative tolerance for all numerical comparisons, By
        default the value returned by the `RTOL` function is used.

    traceback : bool, optional
        If True then print a traceback highlighting where the two
        instances differ.

:Returns: 

    out : bool
        Whether or not the two instances are equal.

**Examples**

'''
        if self is other:
            return True
        
        # Check that each instance is the same type
        if self.__class__ != other.__class__:
            if traceback:
                print("%s: Different types: %s, %s" %
                      (self.__class__.__name__,
                       self.__class__.__name__,
                       other.__class__.__name__))
            return False
        #--- End: if

#        # Check the attributes
#        attrs = set(self.__dict__)
#        if attrs != set(other.__dict__):
#            if traceback:
#                print("%s: Different attributes: %s" %
#                      (self.__class__.__name__,
#                       attrs.symmetric_difference(other.__dict__)))
#            return False
#        #--- End: if

        if rtol is None:
            rtol = RTOL()
        if atol is None:
            atol = ATOL()

#        for attr in attrs.difference(('_dict',)):
#            x = getattr(self, attr)
#            y = getattr(other, attr)
#            if not equals(x, y, rtol=rtol, atol=atol, traceback=traceback):
#                if traceback:
#                    print("%s: Different '%s' attributes: %s, %s" %
#                          (self.__class__.__name__, attr, x, y))
#                return False
#        #--- End: for

        # Check that the keys are equal
        if set(self) != set(other):
            if traceback:
                print("%s: Different keys: %s" %
                      (self.__class__.__name__,
                       set(self).symmetric_difference(other)))
            return False
        #--- End: if

        # Check that the key values are equal
        for key, value in self.iteritems():
            if not equals(value, other[key], rtol=rtol, atol=atol,
                          traceback=traceback):
                if traceback:
                    print("%s: Different '%s' values: %s, %s" %
                          (self.__class__.__name__, key,
                           repr(value), repr(other[key])))
                return False
        #--- End: for
                
        # Still here?
        return True
    #--- End: def

    def get_keys(self, regex=None):
        '''

Return a list of the key names which match a regular expression.

:Parameters:

    regex : str, optional
        The regular expression with which to identify key names. By
        default all keys names are returned.

:Returns: 

    out : list
        A list of key names.

**Examples**

>>> d.keys()
['dim2', 'dim0', 'dim1', 'aux0', 'cm0']
>>> d.get_keys()
['dim2', 'dim0', 'dim1', 'aux0', 'cm0']
>>> d.get_keys('dim')
['dim2', 'dim0', 'dim1']
>>> d.get_keys('^aux|^dim')
['dim2', 'dim0', 'dim1', 'aux0']
>>> d.get_keys('dim[1-9]')
['dim2', 'dim1']

'''
        if regex is None:
            return self.keys()

        keys = []        
        for key in self:
            if re_search('%s' % regex, key):
                keys.append(key)
        #--- End: for

        return keys
    #--- End: def

    def has_key(self, key):
        '''

Return true if and only if the dictionary contains the given key.

:Parameters:

    key : hashable object
        The key.

:Returns:

    out : bool

**Examples**

>>> d.keys()
['key1', 3, ('a', 1)]
>>> d.has_key(3)
True
>>> d.has_key('key9')
False

'''
        return self._dict.has_key(key)
    #--- End: def

#--- End: class


#class foo(CfDict):
#    '''
#'''
#    def __init__(self, *args, **kwargs):
#        '''
#
#**Initialization**
#
#The `dict_name` attribute is initialized to the class name.
#
#:Parameters:
#
#    args, kwargs
#        Keys and values are initialized exactly as for a built-in
#        dict.
#
#**Examples**
#
#
#'''
#        super(foo, self).__init__(*args, **kwargs)
#
#        self.dict_name =  self.__class__.__name__
#        '''
#The name of the dictionary.
#'''
#   #--- End: def    
#
##--- End: class
#
#
#class prop(foo):
#    '''
#'''
#    pass
##--- End: class
#
#
#class attr(foo):
#    '''
#'''
#    pass
##--- End: class
#
#
#class coord(foo):
#    '''
#'''
#    pass
##--- End: class
#
#
#class cellsize(foo):
#    '''
#'''
#    pass
##--- End: class

def RTOL(*rtol):    
    '''

Return or set the default value of relative tolerance for testing
numerically tolerant equality.

:Parameters:

    rtol : int, optional
        The new value of relative tolerance.

:Returns:

    out : float or None    
        If `rtol` was not set return the existing value of relative
        tolerance, otherwise return `None`.

**Examples**

>>> cf.RTOL()
1.0000000000000001e-05
>>> cf.RTOL(1e-10)
>>> cf.RTOL()
1e-10

'''
    if rtol:
        CONSTANTS['RTOL'] = rtol[0]
    else:
        return CONSTANTS['RTOL']
#--- End: def

def equals(x, y, rtol=None, atol=None, traceback=False):
    '''

True if two objects are logically equal, False otherwise.

If the first argument, `x`, has an `equals` method then it is used,
and in this case ``equals(x, y)`` is equivalent to ``x.equals(y)``.

:Parameters:

    x, y :
        The objects to compare for equality.

    atol : float, optional
        The absolute tolerance for all numerical comparisons, By
        default the value returned by the `ATOL` function is used.

    rtol : float, optional
        The relative tolerance for all numerical comparisons, By
        default the value returned by the `RTOL` function is used.

    traceback : bool, optional
        If True then print a traceback highlighting where the two
        objects differ.

:Returns: 

    out : bool
        Whether or not the two objects are equal.

**Examples**

>>> x
<CF Field: rain(10,20)>
>>> cf.equals(x,x)
True

>>> cf.equals(1.0, 1.0)
True
>>> cf.equals(1.0, 33)
False

>>> cf.equals('a', 'a')
True
>>> cf.equals('a', 'b')
False

>>> type(x), x.dtype
(<type 'numpy.ndarray'>, dtype('int64'))
>>> y=x.copy()
>>> cf.equals(x, y)
True
>>> cf.equals(x, x+1)
False

>>> class A(object): pass
>>> a=A()
>>> b=A()
>>> cf.equals(a, a)
True
>>> cf.equals(a, b)
False

'''

    if rtol is None:
        rtol = RTOL()
    if atol is None:
        atol = ATOL()

    if hasattr(x, 'equals') and callable(x.equals):
        # x has a callable equals method
        return x.equals(y, rtol=rtol, atol=atol, traceback=traceback)

    else:
        if hasattr(x, '__iter__') or hasattr(y, '__iter__'):
            # x or y is a sequence  
            if not isinstance(x, numpy.ndarray):
                x = numpy.ma.asarray(x)
            if not isinstance(y, numpy.ndarray):
                y = numpy.ma.asarray(y)
                   
            if x.shape != y.shape:
                return False

            return array_allclose(x, y, rtol=rtol, atol=atol)
        else:
            try:
                # x and y are both numbers
                return abs(x-y) <= atol + rtol*abs(y)
            except TypeError:
                # At least one of x and y is not a number
                return x == y
        #--- End: if
    #--- End: if

#--- End: def

def dump(x, **kwargs):
    '''
    
Print a description of an object to stdout.

If the object has a `dump` method then this is used to create the
output. In this case the arguments are passed to the `dump`
method. Otherwise the arguments are ignored and ``str(x)`` is printed.

:Parameters:

    x : object, optional
        The object to print.

    kwargs :

:Returns:

    None

**Examples**

>>> cf.dump(x)
>>> cf.dump(f, id='field2')

'''
    if hasattr(x, 'dump') and callable(x.dump):
        print x.dump(**kwargs)
    else:
        print x
#--- End: def

def array_allclose(a, b, rtol=None, atol=None):
    '''

True if two arrays have the same shape and elements to within
numerical tolerance, False otherwise.

The tolerance values are positive, typically very small numbers. The
relative difference (`rtol` * abs(`b`)) and the absolute difference
`atol` are added together to compare against the absolute difference
between `a` and `b`.

:Parameters:

    a, b : array_like
        Input arrays to compare.

    atol : float, optional
        The absolute tolerance for all numerical comparisons, By
        default the value returned by the `ATOL` function is used.

    rtol : float, optional
        The relative tolerance for all numerical comparisons, By
        default the value returned by the `RTOL` function is used.

:Returns:

    out : bool
        Returns True if the arrays are equal.

**Examples**

>>> array_allclose([1, 2], [1, 2])
True
>>> array_allclose(np.array([1, 2]), np.array([1, 2]))
True
>>> array_allclose([1, 2], [1, 2, 3])
False
>>> array_allclose([1, 2], [1, 4])
False

'''      
    try:
        return numpy.ma.allclose(a, b, rtol=rtol, atol=atol)
    except (IndexError, NotImplementedError):
        return numpy.ma.all(a == b)
#--- End: def

def _files(openfile):
    _open_files.append(openfile)
#--- End: def

def close():
    for f in _open_files:
        f.close()
    _open_files[:] = []
#--- End: def

def parse_indices(data, indices):
    '''

'''
    parsed_indices = []

    if not isinstance(indices, tuple):
        indices = (indices,)

    # Initialize the list of parsed indices as the input indices with any
    # Ellipsis objects expanded
    length, dims = len(indices), len(data.shape)
    for index in indices:
        if index is Ellipsis:
            parsed_indices.extend([slice(None)] * (dims-length+1))
            length = len(parsed_indices)
        else:
            parsed_indices.append(index)
    #--- End: for

    ndim = data.ndim

    if ndim and len(parsed_indices) > ndim:
        raise IndexError("Invalid indices %s for array with shape %s" %
                         (parsed_indices, data.shape))

    if len(parsed_indices) < ndim:
        parsed_indices.extend([slice(None)]*(ndim-len(parsed_indices)))

    if not ndim and parsed_indices:
        # If data is scalar then allow it to be indexed with an
        # equivalent to [0]
        if (len(parsed_indices) == 1 and
            parsed_indices[0] in (0, -1, slice(0,1), 
                                  slice(-1,None,-1), slice(None, None, None))):
            parsed_indices = []
        else:            
            raise IndexError(
"Scalar array can only be indexed with (), Ellipsis or an equivalent to 0")

    #--- End: if

    for i, (index, size) in enumerate(zip(parsed_indices, data.shape)):

        if isinstance(index, slice): 
            start, stop, step = index.indices(size)
            if (start == stop or
                (start < stop and step < 0) or
                (start > stop and step > 0)):
                raise IndexError("Invalid indices %s for array with shape %s" %
                                 (parsed_indices, data.shape))
            if step < 0 and stop < 0:
                stop = None
            index = slice(start, stop, step)            
        elif isinstance(index, (int, long)):
            if index < 0: 
                index += size
            index = slice(index, index+1, 1)
        else:    
            if getattr(getattr(index, 'dtype', None), 'kind', None) == 'b':
                # Convert booleans to +ve integers
                index = list(numpy.where(index)[0])
            else:
                # Convert negative integers to non-negative integers
                index = [(x+size if x<0 else x) for x in index]
        
            if len(index) == 1:
                # Convert a single element list to a slice object
                index = index[0]
                index = slice(index, index+1, 1)
            else:                    
                # Try to find a slice object equivalent to the list
                step = index[1] - index[0]
                if step > 0:
                    start, stop = index[0], index[-1]+1
                elif step < 0:
                    start, stop = index[0], index[-1]-1

                if index == range(start, stop, step):
                    # Replace the list with a slice object
                    if stop < 0:
                        stop = None
                    index = slice(start, stop, step)
            #--- End: if
        #--- End: if
                    
        parsed_indices[i] = index    
    #--- End: for
    
    return parsed_indices
#--- End: def

def subspace_array(array, indices):
    '''

Subset the input numpy array with the given indices. Indexing is similar to
that of a numpy array. The differences to numpy array indexing are:

1. An integer index i takes the i-th element but does not reduce the rank of
   the output array by one.

2. When more than one dimension's slice is a 1-d boolean array or 1-d sequence
   of integers then these indices work independently along each dimension
   (similar to the way vector subscripts work in Fortran).

indices must contain an index for each dimension of the input array.
'''
    gg = []
    for i, x in enumerate(indices):
        if not isinstance(x, slice):
            gg.append(i)
    #--- End: for

    len_gg = len(gg)
    if len_gg > 1:
        # Slice the dimensions one at a time
        indices = list(indices)
        for axis in gg:
            array = numpy.ma.take(array, indices[axis], axis=axis)
            indices[axis] = slice(None)
        if len_gg < len(indices):
            array = array[tuple(indices)]
    else:
        # Slice all dimensions at the same time
        array = array[tuple(indices)]

    return array
#--- End: def

def iterindices(location):
    '''

Return an iterator over indices ...

:Parameters:

    location: sequence

:Returns:

    out : generator
        An iterator over the element's indices

**Examples**

>>> for index in iterindices(([1,3], [0, 1], [3, 6])):
...     print index
...
(1, 0, 3)
(1, 0, 4)
(1, 0, 5)
(2, 0, 3)
(2, 0, 4)
(2, 0, 5)

>>> for index in iterindices([]):
...     print index
...
()

>>> for index in iterindices([(0, n) for n in [2, 2]])
...     print index
...
(0, 0)
(0, 1)
(1, 0)
(1, 1)

'''
    indices = [xrange(*r) for r in location]
    for index in itertools_product(*indices):
        yield index
#--- End: def

def FM_THRESHOLD(*new_minncfm):
    '''

Return or set the minimum amount of memory to be kept free as a
temporary work space.

The amount is returned as a number of kibibytes of memory, but set as
a number of chunks.

:Parameters:

    new_minncfm : int, optional
        The number of chunks to be kept free as a temporary work space.

:Returns:

    out : float or None    
        If new_minncfm was not set return the existing temporary work
        space size in kibibytes, otherwise return `None`.

**Examples**

>>> cf.FM_THRESHOLD()
1024000.0
>>> cf.FM_THRESHOLD(20)
>>> cf.FM_THRESHOLD()
2048000.0

'''
    if not new_minncfm:
        return CONSTANTS['FM_THRESHOLD']
    
    minncfm = new_minncfm[0]
    CONSTANTS['MINNCFM']      = minncfm
    CONSTANTS['FM_THRESHOLD'] = (CONSTANTS['CHUNKSIZE']/1024.0) * minncfm
#--- End: def

def CHUNKSIZE(*new_chunksize):
    '''

Return or set the chunk size for data storage and processing.

When setting the chunk size, the amount of minimum amount of memory to
be kept free as a temporary work space is also updated.

:Parameters:

    new_chunksize : int, optional
        The new chunk size in bytes.

:Returns:

    out : int or None   
        If new_chunksize was not set then return the existing chunk
        size in bytes, otherwise return `None`.

**Examples**

>>> cf.CHUNKSIZE()
104857600
>>> cf.CHUNKSIZE(2**30)
>>> cf.CHUNKSIZE()
1073741824

'''
    if not new_chunksize:
        return CONSTANTS['CHUNKSIZE']
    
    CONSTANTS['CHUNKSIZE'] = new_chunksize[0]

    FM_THRESHOLD(CONSTANTS['MINNCFM'])
#--- End: def

# --------------------------------------------------------------------
# Set the default values of MINNCFM, CHUNKSIZE and FM_THRESHOLD
# --------------------------------------------------------------------
CONSTANTS['MINNCFM'] = 10
CHUNKSIZE(104857600)       #  = 100 mebibtyes

